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A two- component Bose-Einstein condensate rotating in a toroidal trap is investigated. The topo- 
logical constraint depends on the density distribution of each component along the circumference of 
the torus, and therefore the quantization condition on the circulation can be controlled by changing 
the miscibility using the Feshbach resonance. We find that the system exhibits a variety of dynamics 
depending on the initial angular momentum when the miscibility is changed. 
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I. INTRODUCTION 

A Bose-Einstein condensate (BEC) in a toroidal 
trap [1-4] is an ideal system to study fascinating proper- 
ties of a superfluid, such as persistent flow symmetry 
breaking localization [E|-[8|, and various rotating states 
arising from quantized circulations 0, [To| • 

The quantization of circulation in superfluids origi- 
nates from the single-valuedness of the wave function, 
i.e., the change in the phase of the wave function must 
be an integer multiple of 2tt along a closed path [ll|, [l2| • 
An angular momentum of a BEC in a toroidal trap is 
therefore quantized if the density is uniform. However, if 
a density vanishes at a part of the circumference of the 
torus, the phase can jump at the density defect Il3j |, 
and the system is allowed to rotate with an arbitrary 
circulation. 

In the present paper, we consider a system of a two- 
component BEC rotating in a toroidal trap [3, [H| . Let 
us consider a situation in which the repulsive interaction 
separates the two components along the circumference of 
the torus as illustrated in Fig. [TJ Since the closed path 
along the torus for one component is blocked by the other 
component, the quantization condition is not imposed on 
the circulation around the torus, and the system can ro- 
tate at an arbitrary frequency. If we decrease the inter- 
component repulsion from this initial state, the two com- 
ponents mix, and the density distribution of each com- 
ponent becomes multiply connected, and consequently, 
the quantization condition is imposed on the circulation 
along the torus for each component. During this dy- 
namics, the total angular momentum must be conserved. 
Nontrivial dynamics is thus expected due to the interplay 
between the change of the quantization condition on the 
circulation and the angular momentum conservation. 

This paper is organized as follows. Section [II] stud- 
ies one-dimensional (ID) ring geometry. Section III Al 
gives formulation of the problem. Section III Bl numer- 
ically investigates the stationary states as a function of 
the intercomponent repulsion for a given angular momen- 
tum, and Sec. Ill CI calculates the Bogoliubov spectrum. 
Section III PI compares analytic results with numerical 
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FIG. 1: (color online) Schematic illustration of a phase- 
separated two-component condensate in a toroidal trap. 



ones. Section IIIII demonstrates the dynamics of a two- 
component Rb BEC in a 3D toroidal trap. Section llVl 
gives conclusions to this study. 



II. ONE DIMENSIONAL RING 

A. Formulation of the problem 

In the mean-field approximation, a two-component 
BEC in a frame rotating at a frequency ft around the 
z axis is described by the Gross-Pitaevskii (GP) equa- 
tions, 
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where ipj is the macroscopic wave function for the jth 
component (j — 1, 2) normalized as J \ipj\ 2 dr — Nj with 
Nj being the number of atoms, rrij is the atomic mass, Vj 
is the potential, and L z = —ifid/dO is the z component 
of the angular momentum operator with = dLig(x + iy). 
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The interaction parameters are given by 



g jr = 27rh 2 a jr [m j 1 + m/) , 



(2) 



where a^y is the s- wave scattering length between atoms 
in components j and j' . The condition for the phase 
separation in an infinite system is given by gii#22 < 9\2 
for g jr > p3|. 

In this section, for simplicity, we reduce the problem to 
ID with a periodic boundary condition. We also assume 
rrii = m 2 = m, V\ = V2, and gn = #22- Assuming that 
the system is tightly confined in the torus and neglecting 
excitations in the radial and z directions, we obtain the 
normalized GP equation as (j ^ j f ) 
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where the wave function ipj(0) is normalized as 



/hM 2 d0 = N j /(N 1 +N 2 ) 



r M/{mR 2 ), Q 



mR 2 Q/h, and g^y = mR 2 gjj> p av /ft 2 with R being the 
radius of the torus and p av the averaged density. 



B. Stationary state with fixed angular momentum 
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We consider a situation in which a phase separated 
stationary state in a rotating frame is prepared and then 
the inter component interaction #12 > is adiabatically 
decreased, during which the angular momentum of the 
system is conserved. We therefore seek the stationary 
state of Eq. (j3j) for a given angular momentum as a func- 
tion of #12. 

We use the imaginary time propagation method, in 
which i on the left-hand side of Eq. (j3j) is replaced by 
— 1. We add —jlj\ipj\ 2 to the right-hand side of Eq. (|3j) 
and control jlj and Q in the imaginary time propagation 
such that the norm J \ipj\ 2 d0 = nj in each component 
and the total angular momentum, 




are kept constant. The wave function thus converges to 
the stationary state with given norms nj and an angular 
momentum L. 

Figure [2] (a) shows the density and phase profiles of the 
stationary states of Eq. (|3j), where the populations are 
77-1 =0.9 and ri2 = 0.1 and the total angular momentum 
is fixed to L = 0.3. In Fig. [2j we start from §12 = 1200 
and slowly decrease §12 in the imaginary time propaga- 
tion. At §12 = 1200 [top panels of Fig. [2] (a)], the two 
components are phase separated. Component 2 breaks 
component 1 at 6 = 0, where the phase of component 1 
changes rapidly. As g±2 is decreased, the two components 
are mixed, and the density distribution of each compo- 
nent becomes uniform for §12 < 980 [bottom panels of 



FIG. 2: (color online) (a) Density distributions (left panels) 
and phase profiles (right panels) of the stationary states of 
Eq. © for 012 = 1200, 1000, and 900. (b) Normalized angular 
momenta L\ (solid curve), L2 (dashed curve), and rotation 
frequency (l (circles) of the stationary states as functions of 
g\2- In (a) and (b), gn = gii — 1000 and the total angular 
momentum is fixed to L — 0.3. The populations are m =0.9 
and 712 = 0.1. 

Fig. |5](a)]. We note that the circulation, 

changes due to the constraint of L = 0.3 being fixed. As 
§12 is decreased, component 1 goes to the uniform state 
with Ti = 0, and in order to maintain L = 0.3 the circu- 
lation of component 2 increases to T2 = 3, which satisfies 
n\Yi + 712^2 = L. Thus L 2 increases with a decrease in 
L\ as shown by the solid and dashed curves in Fig.[2](b). 
The plots in Fig. [2] (b) show the rotation frequency ft of 
the frame in which the density distributions are station- 
ary. We find that Q coincides with Z/2/^2, indicating the 
rigid-body rotation of component 2 for §12 < 980. 

Figure [3] shows the case of L = 0.5. Unlike in Fig.[2j the 
dip in |^i | 2 persists until g±2 = 900, which is occupied by 
the localized component 2. The stability of this structure 
is due to the fact that arg^i jumps by 7r at = and the 
stable dark soliton is formed in component 1. As g±2 is 
decreased, the tail of |^2| 2 spreads and the quantization 
condition is imposed on the circulation [bottom panels of 
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FIG. 3: (color online) (a) Density distributions (left panels) 
and phase profiles (right panels) of the stationary states of 
Eq. J3J. (b) Normalized angular momenta L± (solid curve), 
L2 (dashed curve), and rotation frequency Q (circles) of the 
stationary states as functions of gi2. The total angular mo- 
mentum is fixed to L = 0.5. The other conditions are the 
same as those in Fig. [2] 



FIG. 4: (color online) (a) Density distributions (left panels) 
and phase profiles (right panels) of the stationary states of 
Eq. (|3]). (b) Normalized angular momenta L\ (solid curve), 
L2 (dashed curve), and rotation frequency Q (circles) of the 
stationary states as functions of g\2- The total angular mo- 
mentum is fixed to L — 0.7. The other conditions are the 
same as those in Fig. [2] 



Fig. [3] (a)], which is the reason for the change in Lj at 
012 ^950 in Fig. [3](b). 

Figure |4] shows the case of L = 0.7. As g\2 is decreased, 
ipi goes to a uniform state with circulation Ti = 1 [bot- 
tom panels of Fig. [4] (a)], which has an angular momen- 
tum L\ = n\ = 0.9. As a consequence, L2 must be 
L — Li = —0.2 and therefore Y2 becomes — 0.2/n2 = —2. 
Because of L2 — n2^, the rotation frequency of the sys- 
tem decreases with gi2. This indicates that the rotation 
of the system slows down, stops, and counterrotates in 
real time propagation as gi2 is decreased adiabatically. 
This interesting behavior is due to the interplay between 
the quantization of circulation and the angular momen- 
tum conservation. 

We note that the density distributions in Fig. [2] are 
exactly the same as those in Fig. 2J This is because the 
GP equation (j3j) is invariant under the transformation, 

^ e m tf, (6a) 
£-tt, (6b) 

where i is an integer. By this transformation, the angular 



momentum is changed as 

Lj^inj-Lj. (7) 

The results in Fig.|4]agree with those obtained from Fig. [2] 
by the transformation (|6|) and © with 1 = 1. 

C. Bogoliubov analysis 

If the time scale of the change in g\2 is much longer 
than the inverse of the lowest excitation frequency, the 
wave function follows the stationary states shown in 
Figs. [2][31 In order to find the adiabatic condition, we 
perform the Bogoliubov analysis. 

We separate the wave function as 

^(0,T) = [^ j (9)+5^ j (0,T)]e- i ^, (8) 

where ^ j is a stationary state of Eq. ([3]) and <pj is a small 
deviation. Substituting Eq. (j5J) with 

5^(0, r) = Uj{e)e-^ T +v^{0)e^ T (9) 
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D. Analysis for rotation frequency 
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FIG. 5: (color online) The lowest excitation frequency cDo for 
the stationary state with fixed L as a function of g\2. For 
L — 0.3 and L — 0.7 (solid curve), the excitation frequencies 
are exactly the same. The parameters are gn = #22 = 1000, 
n\ — 0.9, and 712 = 0.1. After the transition to the uniform 
state in the solid curve, Q is fixed. 



into Eq. (j3]) and neglecting the second and higher orders 
of Uj and Vj , we obtain the Bogoliubov-de Gennes equa- 



tion (j^f), 
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We numerically diagonalize Eq. (fTD|) . 

Figure [5] shows the lowest excitation frequency uoq for 
L = 0.3, 0.5, and 0.7. The excitation frequencies for 
L = 0.3 and 0.7 are exactly the same, since both states 
are connected to each other by Eqs. ([6]) and ([7j). We find 
that ujq decreases with §12 and reaches (Do = 0, at which 
the adiabatic condition breaks down. For L = 0.3 and 
0.7 (solid curve), this point (gi2 — 980) corresponds to 
the transition to the uniform state, which indicates that 
the localized state cannot transform to the uniform state 
adiabatically. For L = 0.5 (dashed curve), ujo = at 
012 ^ 940. 

For example, let us consider a 87 Rb BEC with density 
10 14 cm -3 in a ring with radius 10 jam. In this case the 
unit of time is u t = mR 2 /h ~ 0.14 s. If the typical time 
scale of change in gi2 is 1 s ~ Wut, the energy gap cDo 
must be > 1 for the adiabatic condition, and hence g\2 
can be decreased to ~ 990 (~ 960) for the solid (dashed) 
curve in Fig. [5] 



In Sec. IIIB| we showed that the rotation frequency ft 
changes as the intercomponent interaction gi2 is changed. 
In this subsection, we derive an analytic expression of the 
rotation frequency 

We write the wave function as 



M0) = \^P 3 {0)e 



ut>j(e) 



(11) 



where pj{6) and (j)j{6) are real functions. We assume 
that p2(0) is localized near = and define #0 at which 
p2 decays to zero p2(±#o) — 0. We also assume that the 
total density is almost constant, pi{0) J rp2{0) ~ (27r) _1 = 
po, and hence pi(0) falls only around = and pi(0) ~ 
po for \0\ > 0q (— 7r < < 7r). We write the phases as 



h{9) = n e + f(e), 
fo(9) = fie, 



(12a) 
(12b) 



where Qq is a constant determined later. The function 
f(0) changes only around = and satisfies f'{0) ~ 
for \0\ > 0q . From the single- valuedness of the wave 
function, Eq. (|12a|) gives 



2tt^ 



f{0)dO = 2£tt, 
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where I is an integer. Since p2{0) vanishes for \0\ > #o 5 
there is no constraint on Eq. (|12bj) for \0\ > 0q. Using 
dpi/dt = and Eq. ©, we find 



Ji(0) - Qpi(0) = const., 
where Jj is the flux given by 

From Eqs. (fT2b|) . (JUJ), and (JUJ), we obtain 

Po 



/'(6») = -(fi-fi ) 
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The angular momentum of each component in Eq. (j4| is 
rewritten as 



/7T f'TT 
p 1 (9)</>' 1 (9)d9 = / pi( 
-7T J —IT 



(9) n + f(9) 



L 2 = 7l2&, 
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where we used Eqs. ([T2]) and ([16]). Equation (|17b|) is 
consistent with the agreement between the solid curves 
and circles in Figs[2]|4j Equation ([T7]) gives 



L = L\ + L'2 = ho- 
using Eqs. $13]), ([16]), and ([11]), we obtain 



n = i + (l ■ 
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FIG. 6: (color online) Relation between the peak 
density of component 2 and the rotation frequency for 
(L,n 2 ) = (0.3,0.1) (circles), (0.1,0.1) (triangles), and 
(0.4,0.2) (squares). The solid curve shows Eq. ([2T]) . 
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We assume Gaussian density distributions as 

Pl (0) = p - P r k e- e2/a \ (20a) 
P2{9) = p p 2 eak e- e2/a \ (20b) 

where p^ eak = ^/(v^) is the peak density of compo- 
nent 2. The Gaussian is assumed to be very narrow, 
<j <C 1, so that the wave functions smoothly connect at 
= ±7r. Substituting Eq. (j20j) into Eq. ([T9]) . we obtain 



ft = * + (L - i) 



1 



n 2 -^gi/2 

P2 



/ peak \ 



where the function is defined by 

OO 

9i/2(z) = Y.TU2 



k=l 



(21) 
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Figure [6] shows n2^l/ L — 1) as a function of p^ ea / po = 
27rp2 eak - The plots show the numerical results obtained 
by solving Eq. ((3]) and the solid curve shows Eq. (|2T|) 
with t = 0, i.e., 
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We find that the plots with different parameters fit the 
universal curve given by Eq. ([23]) . According to Eqs. (|6j) 
and ([7]) , numerical results with L — ^ i — L also agree well 
with Eq. (jUJ). 



III. THREE DIMENSIONAL TOROIDAL TRAP 

We perform a numerical calculation of full 3D real- 
time propagation in a realistic situation. We consider 
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FIG. 7: (color online) Column density profiles f\ipj\ 2 dz for 
(a) (a, t) = (0.3, 0), (b) (0.3, 60 ms), and (c) (0.5, 80 ms). The 
field of view is 25 x 25 fim. Density profiles of component 2 
at r± = R and z — for (d) a — 0.3, (e) a — 0.5, and (f) 
a — 0.7. The density is normalized by N/R 2 in (a)-(c) and 
by N/R 3 in (d)-(f). 



a BEC of 87 Rb atoms in the hyperfine states \F, mp) = 
1 2, 2) and |2, 1), which we refer to as components 1 and 2, 
respectively. The scattering lengths between component 
1 is an = CL4 and that between component 2 is CL22 = 
(3<22 + 4<24)/7, where as is the s-wave scattering length 
between two atoms with total spin S. Their difference is 
&11 — &22 = 3(<24 — d2)/7 ~ 2.98<2b with <2b being the Bohr 
radius, where we used the value measured in Ref. [l7| . 
We therefore use the scattering lengths an = lOOae and 
&22 = (100 — 2.98)ae. The intercomponent scattering 
length ai2 is assumed to be variable. The number of 
atoms is TV = N ± + N 2 = 2 x 10 5 with N ± /N 2 = 9 - We 
employ a toroidal-shaped trap as 



V X =V 2 



1 



-muj 



(r ± - Rf + z 2 



(24) 



with a frequency uj = 2ir x 1 kHz and a radius R — 10 
/im, where r± = (x 2 + y 2 ) 1 ^ 2 . 
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We first prepare the nonrotating ground state ^j(r) 
for a\2 = lOOae by the imaginary-time propagation of 
Eq. (pQ), which is phase separated as shown in Fig. [71 (a). 
We then imprint the phase as ^ j(r)e %ote with < 6 < 2tt 
for j = 1 and — ir < < ir for j = 2, where a real param- 
eter a determines the initial angular momentum. In an 
experiment, such phase imprinting is done by, e.g., spa- 
tially modulated laser fields [18|. If a is not an integer, 
the phase of component 1 (2) jumps at 6 = (0 = 7r), 
at which the density vanishes. Using this initial state, 
we study the dynamics of the system by solving Eq. (pQ) 
in a laboratory frame (Q = 0) with the pseudospectral 
method. The intercomponent scattering length is de- 
creased as 

a 12 (t) = 100(1 - 0mt/t d )a B (25) 

with td = 1 s. 

Figures [7] (d)-[7] (f) show the dynamics of the density 
profile of component 2 on the circumference of the torus, 
\ip 2 (r± = R,z = 0,t)| 2 . For t < 300 ms, the rotation 
frequencies are almost constant. At a later time, the 
rotation of the system accelerates for a = 0.3 [Fig. [3(d)], 
and reverses for a — 0.7 [Fig. [71 (f)]. For a = 0.5, the 
solitonic structure remains until t = 800 ms as shown 
in Fig. [7] (c). These behaviors are consistent with the 
ID results in Figs. EJU Thus, we have shown that the 
rotation properties studied for a ID ring in Sec. [II] can 
be observed in a 3D toroidal geometry experimentally. 

IV. CONCLUSIONS 

We have studied a two-component BEC rotating in a 
toroidal trap, in which the intercomponent interaction 



is controlled. When the two components are phase sepa- 
rated along the circumference of the torus, the system can 
rotate without topological constraint. As the intercom- 
ponent repulsion is decreased and the two components 
mix, the topological constraint is imposed on the phase 
of the wave function. Thus, the interplay between the 
quantization of circulation and the angular momentum 
conservation exhibits nontrivial phenomena. 

We found that as the two components become misci- 
ble, the system goes to the states with circulations Tj de- 
pending on the initial angular momentum L. For a small 
initial angular momentum (e.g., L = 0.3NH), the major 
component 1 goes to Ti =0 while the minor component 
2 goes to T2 > 0, resulting in the acceleration of rotation 
(Fig. [2]). For a larger angular momentum (L = 0.7 Nh), 
the circulations go to Ti = 1 and T2 < 0, and the system 
counterrotates (Fig. [4]). For L = the stable dark 

soliton is generated (Fig. [3]). The Bogoliubov analysis 
gives the adiabatic condition for the change in the inter- 
component repulsion (Fig. [5]), and the Gaussian analysis 
gives the expression of the rotation frequency (Eq. (|2T|) 
and Fig. [6]). The full 3D numerical analysis has shown 
that the predicted phenomena can be observed in a real- 
istic situation in an experiment. 
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